Finite thermal conductivity in Id lattices 
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We discuss the thermal conductivity of a chain of coupled rotators, showing that it is the first 
example of a Id nonlinear lattice exhibiting normal transport properties in the absence of an on-site 
potential. Numerical estimates obtained by simulating a chain in contact with two thermal baths at 
different temperatures are found to be consistent with those ones based on linear response theory. 
The dynamics of the Fourier modes provides direct evidence of energy diffusion. The finiteness of 
the conductivity is traced back to the occurrence of phase-jumps. Our conclusions are confirmed by 
the analysis of two variants of this model. 



PACS numbers: 

The understanding of heat conduction in insulating 
solids is a long-standing problem even in the simple con- 
text of Id systems. Rigorous studies have shown that the 
thermal conductivity k diverges (in the thermodynamic 
limit) in homogeneous ||l[ and disordered |g] harmonic 
chains. Accordingly, it has become transparent that the 
nonlinear character of the microscopic dynamics is a nec- 
essary ingredient for the emergence of the Fourier heat 
conduction law at the macroscopic level. 

The typical Hamiltonians chosen to study this problem 
belong to the general class 
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where rrii represents the mass of the i-th particle , V 
is the potential energy of internal forces and U is an on- 
site potential representing interactions with a possible 
external substrate. 

The first convincing numerical evidence of a finite con- 
ductivity has been provided by the preliminary study of 
the ding-a-ling model M and the successive detailed anal- 
ysis of the modified ding-dong model [p|. Such results 
have supported the conjecture that a strong chaotic dy- 
namics is the crucial requisite for the observation of nor- 
mal transport properties in Id chains. A fortiori, this 
has been assumed to hold in higher-dimensions, where 
scattering effects, due to nonlinear interactions among 
normal modes, should make the diffusion process, at the 
basis of standard transport phenomena, even more effi- 
cient than in Id systems. 

However, the recent study of the more realistic FPU 
/3-model {V{x) = x'^ /2 + a;4/4, [/ = 0) has revealed 
a power-law divergence of the thermal conductivity as 
K ~ N'^/^ je) (where N is the chain length 0). Such 
an anamalous behaviour has been explained by invok- 
ing the self-consistent mode coupling theory H in the 



description of the effective evolution of long-wavelength 
modes [^. The generality of this result has been con- 
firmed by the study of several models with no on-site 
potential, including the diatomic Toda model [Q that 
was previously believed to yield a finite conductivity. On 
the other hand, numerical results have accumulated in- 
dicating that whenever an external force is added (as in 
the above mentioned ding-a-ling and ding-dong models), 
the conductivity is finite ||ll[. As a consequence, it has 
been suggested that momentum conservation is the key 
ingredient responsible for an anomalous transport even 
in the presence of a chaotic evolution. 

Such a conjecture has been very recently put on more 
firm grounds by Prosen and Campbell |l2|, who have 
shown that anomalous transport necessarily arises when- 
ever two conditions arc fulfilled: (i) the momentum is 
conserved; {ii) the pressure does not vanish in the ther- 
modynamic limit. On the other hand, systems with sym- 
metric potentials, such as the FPU-/3 model, do not be- 
long to this class, since the pressure vanishes. In practice, 
as all momentum-conserving models studied in the liter- 
ature have revealed anomalous transport properties, the 
possibility that the condition on the pressure is redun- 
dant remains open. 

In this Letter, we show that this is not the case by 
analysing some Id models where momentum is conserved 
and yet the conductivity is finite. 

A. The rotator model The simplest example of a 
classical-spin Id model with nearest neighbour interac- 
tions hes in the class (||) with V{x) = 1 — cos(a;) and 
U — Q. This model can be read also as a chain of N 
coupled pendula, where the pi's and the g^'s represent 
action-angle variables, respectively. It has been exten- 
sively studied |13,M as an example of chaotic dynamical 
system |la that becomes integrable both in the small 
and high energy limits, when it reduces to a harmonic 
chain and free-rotators, respectively. In the two inte- 



grable limits, the relaxation to equilibrium slows down 
very rapidly for most of the observables of thermody- 
namic interest (e.g., the specific heat) jlj]. As as conse- 
quence, the equivalence between ensemble and time aver- 
ages is established over accessible time scales only inside 
a limited interval of the energy density e. Here we shall 
discuss heat conduction for values of the energy density 
corresponding to strongly chaotic behaviour. 

B. Numerical analysis of the thermal conductivity 
The most natural and direct way to determine k con- 
sists in simulating a real experiment, by coupling the left 
and right extrema of the chain with two thermal baths at 
temperatures T^ > Tr, respectively. In our simulations 
we have used Nose- Hoover models of thermostats [ [l6[ , 
both because they can be easily implemented (integrat- 
ing the resulting equations with a standard algorithm) 
and because of the smaller finite-size effects (due to the 
unavoidable contact resistance). 

With this setting, a non-equilibrium stationary state 
sets in characterized by a non- vanishing heat flux J ||l^,g| : 
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where fi = -dV{qi+i - qi)/dqi = sin(gi-|.i - qi) is the 
interaction force and ji is the local flux at site i. The 
total heat flux J has to be averaged over a sufficiently 
long time span to get rid of fluctuations and to ensure 
the convergence to the stationary regime. This can be 
tested by monitoring the average heat flux and looking at 
the scale of its fluctuations. As a result, we have verified 
that 2 • 10^ time units are sufficient to guarantee a few 
percents of fluctuations in the worst cases. 

The thermal conductivity is determined by assuming 
the Fourier law, i.e. from the relation J — kVT, where 
VT denotes the imposed thermal gradient. The simu- 
lations have been performed for T^ — 0.55, Tr — 0.35, 
and chain lengths ranging from A^ = 32 to 1024 with fixed 
boundary conditions. The equations of motion have been 
integrated with a 4th-order Runge-Kutta algorithm and 
a time step At — 0.01. The results, reported in Fig. |l| 
clearly reveal a convergence to a K-value approximately 
equal to 7 (see the circles). The dotted line represents the 
best data fit with the function a + b/N: the agreement 
is very good, showing that finite-size corrections to n are 
of the order 0{1/N), as it should be expected because 
of the thermal contacts. However, more important than 
the numerical value of the conductivity is its finiteness in 
spite of the momentum conservation. 

In order to test independently the correctness of our 
results, we have performed direct microcanonical simula- 
tions, which allow determining the thermal conductivity 
through the Green-Kubo formula |18||: 
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where Cj{t) = N{J{t)J{Q)) is the flux autocorrelation 
function at equilibrium and T is the temperature. A cor- 
rect application of the above formula requires fixing the 
energy density e in such a way that the kinetic temper- 
ature (defined as T = (p^), in agreement with the virial 
theorem) is close to the average value of the temperature 
in the previous simulations. The choice s = 0.5 turns 
out to be reasonable, as it corresponds to T w 0.46. In 
the absence of thermal baths, the equations of motion 
are symplectic, so that we have now preferred to use a 
6th order McLachlan-Atela integration scheme 11^ (with 
periodic boundary conditions). 

The correlation function has been computed by ex- 
ploiting the Wiener-Khinchin theorem, i.e. by anti- 
transforming the Fourier power spectrum. The result 
of the time integration is almost independent of A^ for 
A^ > 128. The gray region in Fig. y corresponds to the 
expected value of k taking into account statistical fluctu- 
ations. There is not only a clear confirmation of a finite 
conductivity, but the numerical value obtained with this 
technique is in close agreement with the direct estimates. 

C. Dynamics in the mode space In order to clarify the 
difference between the dynamics of the present model and 
that of FPU-type systems, we have investigated the evo- 
lution of the low-frequency Fourier modes. In Fig. S we 
have reported the power spectra of some long-wavelength 
modes. For the sake of comparison, the same quantities 
are reported for a diatomic FPU-chain, that is charac- 
terized by an anomalous transport. At variance with the 
FPU model, in the rotators there is no sharp peak (which 
is a signal of an effective propagation of correlations M ) . 
Quite differently, the low-frequency part of the spectrum 
is described very well by a Lorentzian with halfwidth 
7 = Dk^ [D K, 4.3). This represents an independent 
proof that energy diffuses, as one expects whenever the 
Fourier law is established. 

D. Temperature dependence of the thermal conductiv- 
ity The most natural question arising from these results 
concerns the reason for the striking difference with other 
symmetric models such as the FPU-/3 system. As long 
as each (gi+i — qi) remains confined to the same valley 
of the potential, there cannot be any qualitative differ- 
ence with the models previously studied in the litera- 
ture. Jumps thorugh the barrier, however, appear to act 
as localized random kicks that contribute to scattering 
the low-frequency modes and thus to a finite conductiv- 
ity. If this intuition is correct, one should find analogies 
between the temperature dependence of the conductiv- 
ity and the average escape time from the potential well. 
To this aim, we have computed k, for different temper- 
ature values by performing microcanonical simulations 
with various energy densities. From the data reported in 
Fig. H, one can notice a divergence for T — > of the type 
K « exp(a/r) with a ~ 1.2. An even more convincing 
evidence of this behavior is provided by the temperature 
dependence of the average escape time (see the triangles 



in Fig. H) with an exponent a w 2. The latter behaviour 
can be explained by assuming that the jumps are the 
results of activation processes. Accordingly, the proba- 
bility of their occurrence is proportional to exp{—AV/T), 
where AV^ is the barrier height and the Boltzmann con- 
stant is fixed equal to 1 (as implicitely done throughout 
the Letter). Since AV = 2, the whole interpretation is 
consistent. Moreover, in the absence of jumps, the de- 
pendence of the conductivity on the length should be the 
same as in FPU-systems, i.e. k « N'^'^. Therefore, a 
low-frequency mode travelling along the chain should ex- 
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perience a conductivity of the order of N , where N 
is the average separation between jumps. Under the as- 
sumption of a uniform distribution of phase jumps, their 
spatial separation is of the same order of their time sep- 
aration, so that we can expect that k w exp[2AU/(5r)]. 
On the one hand, this heuristic argument explains why 
and how such jumps contribute to a normal transport. 
On the other hand, the numerical disagreement between 
the observed and the expected value of the exponent a 
(1.2 vs. 0.8) indicates that our analysis needs refine- 
ments. In fact, we should, e.g., notice that, in the low- 
energy limit, nonlincarities become negligible, implying 
that deviations from the asymptotic law k w 7V^' ^ should 
become relevant. 

E. Further checks In order to test the conjecture that 
jumps are responsible for a normal heat transport, we 
have decided to investigate some other models. First, we 
have considered a double-well potential V{x) = —x^/2 -\- 
x^ /A (the same as in FPU with a different sign for the 
harmonic term). The results of the direct simulations 
are reported in Fig.|l| (see triangles) for a temperature 
corresponding again to a quarter of the barrier height. 
The finiteness of the conductivity ans its numerical value 
is confirmed by the computation of k, through the Green- 
Kubo formula (see the light-grey shaded region). 

Finally, we have considered an asymmetric version 
of the rotator model, namely V{x) = A — cos(a;) + 
0.4sin(2a;), where A is fixed in such a way that the mini- 
mum of the potential energy is zero, and the temperature 
again corresponds to one quarter of the barrier height. 
In this case too, the conductivity is finite, confirming our 
empirical idea that the jumps are responsible for break- 
ing the coherence of the energy flux and, in turn, for the 
finite conductivity. However, from the point of view of 
Ref. Q, this result is quite unexpected, since, in view 
of the model asymmetry, one expects that the pressure 
(f) — {J2fi)/-^ i^ non-zero. Nonetheless, microcanonical 
simulations show that although the distribution of forces 
is definitely asymmetric, their average value is numeri- 
cally 0. This can be understood by noticing that in view 
of the boundedness of the potential, the system cannot 
sustain any compression. Therefore, no contradiction ex- 
ists with Ref. |l|]. 

In conclusion, in this Letter we have reported about 



the first evidence of normal heat transport in Id systems 
with momentum conservation. Such a behaviour appears 
to be connected with jumps between neighbouring po- 
tential valleys. From the dynamical point of view, it is 
natural to ask what are the peculiar properties of such 
jumps that make them so different from other types of 
nonlinear fluctuations that may occur in single-well type 
of potentials. The only clearly distinctive feature that 
we have found is the "hyperbolic" type of behaviour in 
the vicinity of a maximum of the potential which has to 
be confronted with the typical "ellyptic" character of the 
oscillations around the minima. We hope to be able to 
understand in the future whether this is truly the reason 
for the difference in FPU and rotator systems. 
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kind hospitality during the Workshop on "Complexity 
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FIG. 1. Conductivity k versus chain length N as obtained 
from non-equihbrium molecular dynamics. Circles correspond 
to the rotator model with temperatures Tl = 0.55 and 
Tr — 0.35; triangles correspond to the double well poten- 
tial with Tl = 0.04 and Tr = 0.06. The two hues represent 
the best fit with the function a + b/N. The two shaded regions 
represent the uncertainty about the conductivity on the basis 
of the Green-Kubo formula. 




FIG. 2. Power spectra S{w) in arbitrary units of the 
Fourier modes 1,2,4,8,16,32 and 64 (curves from left to right 
in the high frequency 'region) in a chain of length A*' — 1024 
of rotators with energy density with energy density 0.5. The 
curves result from an average over 1000 independent initial 
conditions. In the inset, the same modes are reported for a 
diatomic FPU chain with masses 1, 2 and energy density 8.8. 
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FIG. 3. Thermal conductivity k (labels on the left axis) 
versus the inverse temperature l/T in the rotator model (open 
circles). Triangles correspond to the average time separation 
between consecutive phase jumps (labels on the right axis) in 
the same system. 



